Сравнение алгоритмов LMS и NLMS для адаптивной фильтрации сигналов

Адаптивная фильтрация является ключевой технологией в цифровой обработке сигналов, находящей применение в системах шумоподавления, эхокомпенсации, идентификации систем и многих других областях. В данной статье мы рассмотрим два классических алгоритма адаптивной фильтрации - LMS (Least Mean Squares) и его модификацию NLMS (Normalized LMS), реализованные на языке Engee с использованием пакета EngeeDSP.
LMS (Least Mean Squares) - это стохастический градиентный алгоритм, который минимизирует среднеквадратичную ошибку между желаемым сигналом и выходом адаптивного фильтра. Основное уравнение обновления весов:

w(n+1) = w(n) + μ·e(n)·x(n)

где:

w(n) - вектор весов фильтра на шаге n
μ - шаг адаптации
e(n) - ошибка (разность между желаемым и фактическим выходом)
x(n) - входной вектор

Алгоритм NLMS

NLMS (Normalized LMS) - это модификация LMS, где шаг адаптации нормализуется по мощности входного сигнала:

w(n+1) = w(n) + (μ / (xᵀ(n)x(n) + ε))·e(n)·x(n)

где ε - малая константа для предотвращения деления на ноль.

Анализ первого примера

Первый код сравнивает LMS и NLMS на сигнале с изменяющейся амплитудой:

x_window = 0.05 * randn(1024*10) .* repeat(0.5 .+ 0.5sin.(2π*(0:1024*10-1)/10240), inner=1)
d_window = filt(FIRFilter([0.5, -0.3, 0.2, 0.1, -0.05]), x_window)

Здесь создается:

Шумовой сигнал с амплитудой, модулированной синусоидой
"Желаемый" сигнал, полученный фильтрацией через FIR-фильтр

Анализ проводится путем вычисления частотной характеристики адаптивного фильтра на каждом кадре и визуализации с помощью surface-plot.

using EngeeDSP, DSP, FFTW, Plots

frameSize = 500
Fs = 8192
x_window = 0.05 * randn(1024*10) .* repeat(0.5 .+ 0.5sin.(2π*(0:1024*10-1)/10240), inner=1)
d_window = filt(FIRFilter([0.5, -0.3, 0.2, 0.1, -0.05]), x_window)

HA = EngeeDSP.LMSFilter(
    Algorithm = "LMS",
    FilterLength = 32,
    StepSize = 1.0,
    LeakageFactor = 1.0,
    InitialValueOfFilterWeights = 0,
    AdaptPort = false, 
    ResetPort = "None",
    OutputFilterWeights = true
)

numFrames = length(d_window) ÷ frameSize
response = zeros(512, numFrames)
t_frame = (0:frameSize:length(d_window)-1)/Fs

i = 1
for start_idx in 1:frameSize:length(x_window)-frameSize+1
    end_idx = start_idx + frameSize - 1
    x_frame = x_window[start_idx:end_idx]
    d_frame = d_window[start_idx:end_idx]
    
    if i == 1
        setup!(HA, x_frame, d_frame)
    end
    
    y, e, w = step!(HA, x_frame, d_frame)
    response[:,i] = 20*log10.(abs.(fft([w; zeros(512-length(w))])))
    i += 1
end

p1 = surface(t_frame, Fs*(0:511)/512, response,
    title="LMS Algorithm Response",
    xlabel="Time (s)", ylabel="Frequency (Hz)", zlabel="Power (dB)",
    camera=(-70, 18), zlims=(-80, 2))

HA = EngeeDSP.LMSFilter(
    Algorithm = "Normalized LMS",
    FilterLength = 32,
    StepSize = 0.3,
    LeakageFactor = 1.0,
    InitialValueOfFilterWeights = 0,
    AdaptPort = false,  
    ResetPort = "None",
    OutputFilterWeights = true
)

response = zeros(512, numFrames)

i = 1
for start_idx in 1:frameSize:length(x_window)-frameSize+1
    end_idx = start_idx + frameSize - 1
    x_frame = x_window[start_idx:end_idx]
    d_frame = d_window[start_idx:end_idx]
    
    if i == 1
        setup!(HA, x_frame, d_frame)
    end
    
    y, e, w = step!(HA, x_frame, d_frame)
    response[:,i] = 20*log10.(abs.(fft([w; zeros(512-length(w))])))
    i += 1
end

p2 = surface(t_frame, Fs*(0:511)/512, response,
    title="NLMS Algorithm Response",
    xlabel="Time (s)", ylabel="Frequency (Hz)", zlabel="Power (dB)",
    camera=(-70, 18), zlims=(-80, 2))

plot(p1, p2, layout=(1,2), size=(1000,400))

LMS алгоритм:
- Использует фиксированный шаг адаптации (μ = 1.0)
- Может быть чувствителен к изменению амплитуды входного сигнала
- Быстрее адаптируется, но может быть менее устойчив
NLMS алгоритм:
- Использует нормализованный шаг адаптации (μ = 0.3)
- Автоматически подстраивается под уровень входного сигнала
- Более стабильная работа при изменении амплитуды

Анализ второго примера

Второй код исследует поведение LMS фильтра на стационарных и нестационарных сигналах:

x = 0.05 * randn(1024*10)
d_stat = filt(FIRFilter([0.5, -0.3, 0.2, 0.1, -0.05]), x)
d_nonstat = filt(FIRFilter([0.8, -0.2, 0.0, 0.1, -0.1]), x)

Здесь:

Стационарный случай: неизменная система (один FIR-фильтр)
Нестационарный случай: система меняется (другой FIR-фильтр)

using EngeeDSP, DSP, FFTW

frameSize = 500
Fs = 44100
x = 0.05 * randn(1024*10)
d_stat = filt(FIRFilter([0.5, -0.3, 0.2, 0.1, -0.05]), x)
d_nonstat = filt(FIRFilter([0.8, -0.2, 0.0, 0.1, -0.1]), x)

HA = EngeeDSP.LMSFilter(
    Algorithm = "LMS",
    FilterLength = 32,
    StepSize = 1.0,
    LeakageFactor = 1.0,
    InitialValueOfFilterWeights = 0,
    AdaptPort = false,
    ResetPort = "None",
    OutputFilterWeights = true
)

t_frame = (0:frameSize:length(x)-1)/Fs
response = zeros(512, length(t_frame))
i = 1
for start_idx in 1:frameSize:length(x)-frameSize+1
    end_idx = start_idx + frameSize - 1
    x_frame = x[start_idx:end_idx]
    d_frame = d_stat[start_idx:end_idx]
    
    if i == 1
        setup!(HA, x_frame, d_frame)
    end
    
    y, e, w = step!(HA, x_frame, d_frame)
    response[:,i] = 20*log10.(abs.(fft([w; zeros(512-length(w))])))
    i += 1
end

p1 = surface(t_frame, Fs*(0:511)/512, response,
    title="Stationary Signal Response",
    xlabel="Time (s)", ylabel="Frequency (Hz)", zlabel="Power (dB)",
    camera=(-70, 18), zlims=(-80, 2))

    response = zeros(512, length(t_frame))
i = 1
for start_idx in 1:frameSize:length(x)-frameSize+1
    end_idx = start_idx + frameSize - 1
    x_frame = x[start_idx:end_idx]
    d_frame = d_nonstat[start_idx:end_idx]
    
    y, e, w = step!(HA, x_frame, d_frame)
    response[:,i] = 20*log10.(abs.(fft([w; zeros(512-length(w))])))
    i += 1
end

p2 = surface(t_frame, Fs*(0:511)/512, response,
    title="Non-Stationary Signal Response",
    xlabel="Time (s)", ylabel="Frequency (Hz)", zlabel="Power (dB)",
    camera=(-70, 18), zlims=(-80, 2))

plot(p1, p2, layout=(1,2), size=(1000,400))

Стационарный сигнал:
- Фильтр успешно сходится к оптимальному решению
- Частотная характеристика стабилизируется со временем
Нестационарный сигнал:
- Фильтр пытается адаптироваться к изменяющейся системе
- Характеристика продолжает изменяться, отражая нестационарность

Вывод

Представленные примеры демонстрируют ключевые особенности работы адаптивных фильтров. Выбор между LMS и NLMS зависит от характеристик сигнала и требований к системе. NLMS обеспечивает более стабильную работу при изменении уровня сигнала, в то время как классический LMS может быть предпочтителен в условиях ограниченных вычислительных ресурсов. Анализ нестационарных систем требует особого подхода и может служить темой для дальнейших исследований.
Также исходя их опытов можно сделать следующие выводы:

Для сигналов с изменяющейся амплитудой предпочтительнее NLMS
При работе с нестационарными системами может потребоваться:
- Уменьшение шага адаптации
- Использование алгоритмов с переменным шагом
- Регулярный сброс фильтра
Визуализация частотной характеристики в динамике - мощный инструмент анализа